Tutorial Part IV — Getting started with animal movement

Introduction

To start with, I would like to express my deep thanks to Cedric Scherer for his help in designing and ameliorating the course scripts.

In this tutorial, we will show you how to get started with (real) animal movement data. To put it simple, animal movement data are just points in time and space, i.e. re-locations of an individual that can be provided as a data.frame with four columns. An identifier for the individual is only needed if you have various animals in one data set, i.e. if you have collared several individuals with the same collar or already appended several data sets:

Identifier Timestamp x_coordinate y_coordinate
Fox_1 2022-01-01 08:10:04 13.47027 52.497802
Fox_1 2022-01-01 08:15:04 13.47015 52.497813

From Course2_SpatialR in our teaching collection, you might remember that such data sets can be imported as simple .txt or .csv files. Once you know in which coordinate reference system (CRS) your coordinates are stemming from, you can assign it to the data.frame, thereby you are creating a geo-referenced data set: a SpatialPointsDataframe-Object (if you work with the old R-package {sp}) or an sf-Object (with package {sf}). Looking at these x- and y-coordinates, you might remember that they look like angular units, e.g. decimal degrees, and the x-coordinate refers to longitude and the y-coordinate to latitude. Hence, we are dealing with EPSG-code 4326. Once you have imported and geo-referenced your data, you can plot and thoroughly check the data.

Checking your data is a mandatory prerequisite before any analysis (see Zuur et al. 2010). With your movement data, you might have outliers in space (e.g. when you test your collars in location A which might be hundreds of kilometers from your actual trapping site), your data might not have been regularly sampled because the animal was hiding and there was no access to the satellite (weak GPS-signal) etc etc. This uneven sampling can affect calculations of speed, for example, and therefore a series of packages have been developed to deal with those issues in movement data, like {move} or {amt}.

The most prominent movement analyses comprise home range estimation, calculations of habitat preference, behavioral analyses (e.g. Hertel et al. 2021) and detection of movement syndromes (personality differences) (Michelangeli et al. 2021). Some of the packages listed in the next chapter are designed for these analyses.

In the following, we will step by step start with loading and exploring movement data of the female red fox (Vulpes vulpes) ‘Q von Stralau’, who had a small home range in Berlin, Germany. She was collared in January 2018 and had a very stable daily routine as shown by 4 (20) min relocation intervals. Her location data are stored in the file with the tag number tag5334_gps.txt. For any details on the data please refer to Kimmig (2021) and Scholz (2020).

We will, however, only use the data of one month in this exercise, as data sets quickly get too big. Please note: the courses Course1_IntroR and Course2_SpatialR of our course script collection
are obligatory for this tutorial.

Useful (web)sites and reading

  • For analysis of telemetry data: packages adehabitatHR, adehabitatLT, move, move2, recurse, momentuHMM, moveHMM, ctmm, amt

Methods papers

Example papers

Getting started

To follow the tutorial, you can either clone or download the repository or you create your own R-project, copy the raw data and type the code chunks into an R-Script. Please refer to the section on using R-projects in Course2_RSpatial of our teaching collection.

If you start with your own R-project, I strongly recommend to use the {d6}-package Cedric Scherer provided. This package automatically sets up the ideal folder structure: https://github.com/EcoDynIZW/d6

In any case, the course folder has the following structure:

#.
#└── d6_teaching_collection         – (root folder)
#    ├─── data                      
#       ├─── data_move              – (contains the file with the GPS locations)
#    ├─── docs
#    ├─── output                    – (contains the cleaned files and processed data)
#    ├─── plots
#    ├─── R                         – (contains the R-script)
#    └ d6_teaching_collection.Rproj – (the RStudio-Project file location)


Necessary packages to install and load

We first have to install the packages and load them before we can use the functions we need.

## package names:
pkgs = c("here", "lubridate", "sf", "sp", "ggplot2", "tmap", "circular", "move", 
         "viridis", "devtools", "plotly","amt", "terra","effects","tidyterra",
         "corrplot","ctmm","sjPlot","jtools") 

# install.packages(pkgs) # only run this line once! for installing the packages!
# update.packages()

Tip of the day: If you have already installed some of the packages above, you can first check which ones are already installed, and save the ones not installed in an object called ‘my_packages’ and only install the missing ones:

my_packages <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
my_packages
## character(0)
if(length(my_packages) > 0) install.packages(my_packages)

Now load the packages:

library(here)          ## for easy directory management
library(lubridate)     ## for date handling
library(sf)            ## for handling simple feature objects
library(terra)         ## for handling raster objects
library(ggplot2)       ## for nice static plots
library(tmap)          ## for nice interactive maps
library(circular)      ## for circular stats + plots
library(move)          ## for `move` objects
library(amt)           ## for `move` objects
library(ctmm)          ## for fitting continuous-time movement models
library(effects)       ## for visualizing model results
#library(sjPlot)        ## easy plotting of model results
library(jtools)        ## easy plotting of model results
library(viridis)       ## for perceptually uniform color palettes
library(corrplot)      ## for correlation plot

Set the working environment

Now that we are working inside an R-project, we can use the easy functionality of the {here} package, which is automatically pointing to your project’s root folder, e.g.:

here::here()

Hence, there is no need to use the function setwd() any more. Note: if it does not work, please close RStudio, go to your Explorer and double-click on the .Rproj file. Then, under ‘files’ (usually lower right panel) double-click on the R folder and open the script.

Load data

The movement data are stored in the data-raw subfolder. Let’s check which files are available.

lf <- list.files(path = here("data","data_move"), full.names = TRUE) 
lf
## [1] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.cpg"
## [2] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.dbf"
## [3] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.prj"
## [4] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.qpj"
## [5] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.shp"
## [6] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.shx"
## [7] ".../d6_teaching_collection/data/data_move/geo-raw"                     
## [8] ".../projects_github/d6_teaching_collection/data/data_move/KettleHoles.txt"             
## [9] ".../d6_teaching_collection/data/data_move/tag5334_gps.txt"

The output lists two results, there is another subfolder, and many files. E.g., the element 7 of the vector lf, lf[7], is a folder. If you already know that your movement data contains e.g. ‘gps’ in its name or is stored as ‘.txt’ or ‘.csv’ files, you can directly search for those files with the pattern argument:

## check the difference, and note: full.names is set to FALSE
thefile <- list.files(path = here("data","data_move"), pattern = "gps", full.names = FALSE)
thefile
## [1] "tag5334_gps.txt"
## ...and here to TRUE
thefullfile <- list.files(path = here("data","data_move"), pattern = "gps", full.names = TRUE)
thefullfile
## [1] "C:/Users/admin/d6_teaching_collection/data/data_move/tag5334_gps.txt""

Now load the data file. Our first fox filename should be tag5334_gps.txt, which we stored under object ‘thefullfilename’

dat_anim <- read.table(file=thefullfile, header=TRUE, fill=TRUE, sep=',')


Data check and cleaning

Checking for missing or incomplete information

Let’s have a look at the data. This is the typical way data are stored on e-obs collars:

dat_anim[1:5,] ## recap: head(dat_anim) also works
##   key.bin.checksum tag.serial.number         start.timestamp longitude latitude
## 1       2135768252              5334 2018-10-30 04:00:00.000  13.47047 52.49499
## 2       2210577241              5334 2018-10-30 04:20:00.000  13.47157 52.49526
## 3       1824790132              5334 2018-10-30 04:40:00.000  13.47103 52.49525
## 4       4079427513              5334 2018-10-30 07:00:00.000  13.47109 52.49526
## 5        380027806              5334 2018-10-30 08:00:00.000  13.47105 52.49527
##   height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 1                   76.3           3      A                   22
## 2                   81.6           3      A                   18
## 3                   66.5           3      A                   53
## 4                   73.5           3      A                   24
## 5                   71.1           3      A                   47
##          timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 1 2018-10-30 04:00:23.000            3604                3389          10
## 2 2018-10-30 04:20:19.000            3616                3333          11
## 3 2018-10-30 04:40:54.000            3624                3374          13
## 4 2018-10-30 07:00:25.000            3659                3491          11
## 5 2018-10-30 08:00:48.000            3660                3470          13
##   speed.over.ground heading.degree speed.accuracy.estimate
## 1              1.64         233.31                    1.67
## 2              4.43         327.02                    3.79
## 3              0.07           0.00                    0.73
## 4              0.01           0.00                    0.65
## 5              0.04           0.00                    0.89
##   horizontal.accuracy.estimate
## 1                         7.94
## 2                        59.65
## 3                        15.36
## 4                         7.17
## 5                         8.45

For now, we will only work with few columns. tag.serial.number refers to the individual identifier (collar ID). There are two columns with timestamps: start.timestamp is the preprogrammed time-interval. Then there is the timestamp.of.fix, which is the real time the GPS-location was recorded. This is usually a bit later, i.e.
( = start.timestamp + used.time.to.get.fix + 1 second), as it takes some time for the collar unit to connect to the satellite. The spatial info is stored in the columns longitude and latitude.

Before you can transform the data.frame dat_anim into a georeferenced spatial object, you need to check whether there are missing locations in your data.frame, otherwise you will get an error message on transformation.

This can happen if no GPS-signal could be recorded. Or - importantly - some collars are only activated when the animal is moving to save battery life. In that case, the missing GPS coordinates would correspond with the last position (= be the same). Depending on which analysis you want to do, e.g. define resting places, you might need to fill the missing positions again.

In our case, there are a lot of missing values in the locations:

## if the latitude-entry is missing, the longitude value will also be missing
## so it is enough to only check the latitude
which(is.na(dat_anim$latitude))
##   [1]   48   65   89   90   91  129  183  185  222  224  226  229  260  266  355
##  [16]  358  415  421  422  423  457  474  510  511  577  601  616  621  626  646
##  [31]  673  674  685  696  721  774  784  796  822  833  842  855  863  907  921
##  [46]  940  993 1006 1038 1070 1077 1078 1088 1126 1150 1159 1183 1188 1215 1238
##  [61] 1268 1272 1304 1330 1480 1516 1550 1552 1559 1600 1637 1680 1717 1721 1736
##  [76] 1738 1742 1772 1820 1824 1830 1867 1884 1973 1997 1999 2020 2023 2024 2025
##  [91] 2026 2064 2075 2076 2109 2162 2234 2313 2329 2374 2379 2396 2419 2453 2455
## [106] 2518 2551 2609 2614 2638 2687 2728 2745 2770 2776 2783 2800 2804 2805 2854
## [121] 2926 2927 2931 2935 2943 2950 2961 2963 2968 2990 2991 3043 3044 3052 3090
## [136] 3102 3114 3172 3225 3226 3305 3328 3361 3372 3373 3386 3409 3498 3516 3517
## [151] 3518 3520 3591 3622 3623 3624

Delete rows with missing spatial info:

dat_anim_na <- dat_anim[!is.na(dat_anim$latitude),] ## alternatively, use complete.cases()

Make the crosscheck if there is missing info in a row in longitude. This should NOT be the case after we had deleted those rows:

which(is.na(dat_anim_na$longitude)) # none
## integer(0)

Checking for coarse spatial outliers

It could happen the collar was tested e.g. in Berlin, but the animal was finally caught and collared far away. Make a quick check whether there are strange locations:

plot(dat_anim_na$latitude ~ dat_anim_na$longitude)

There do not seem to be coarse outliers, data look compact.

Working with the date format

We will now add some additional columns, where separate days (numbered from 1 to 365), the month (from 1 to 12) and the hour of the day are stored (from 0 to 23). Sunset and sunrise can be calculated based on dates and the location (latitude, longitude). This can be done with the package {lubridate}. Check the vignette!

## have a look at the timestamp format, here the first row of the column start.timestamp
dat_anim_na$start.timestamp[1]
## [1] "2018-10-30 04:00:00.000"
## define date-time format - the format is year-month-day_hour:min:sec:
dat_anim_na$start.timestamp <- ymd_hms(dat_anim_na$start.timestamp, tz="Europe/Berlin") 
dat_anim_na$start.timestamp[1]
## [1] "2018-10-30 04:00:00 CET"

Now append the information:

dat_anim_na$yearday <- yday(dat_anim_na$start.timestamp)
dat_anim_na$month   <- month(dat_anim_na$start.timestamp)
dat_anim_na$hour    <- hour(dat_anim_na$start.timestamp)
dat_anim_na$kweek   <- week(dat_anim_na$start.timestamp)
dat_anim_na$date    <- date(dat_anim_na$start.timestamp)
## crosscheck with
# head(dat_anim_na)

In addition, we want to calculate the hours of sunset and sunrise as well as daylength. For this, we need to install a package that is still under development, i.e. which is not on CRAN. We therefore must download and install it locally:

# devtools::install_github("bgctw/solartime") ## run only once to install package from GitHub
library(solartime)

For computing sunset and sunrise, the latitude and longitude must be provided as well:

dat_anim_na$sunrise   <- computeSunriseHour(timestamp = dat_anim_na$start.timestamp,
                                            latDeg = dat_anim_na$latitude,
                                            longDeg = dat_anim_na$longitude)

dat_anim_na$sunset    <- computeSunsetHour(dat_anim_na$start.timestamp,
                                           dat_anim_na$latitude,
                                           dat_anim_na$longitude)

dat_anim_na$daylength <- computeDayLength(dat_anim_na$start.timestamp,dat_anim_na$latitude)

dat_anim_na$daytime   <- computeIsDayByLocation(dat_anim_na$start.timestamp,
                                                dat_anim_na$latitude,
                                                dat_anim_na$longitude)

## check new variables
# head(dat_anim_na)
hist(dat_anim_na$daylength)

unique(dat_anim_na$daytime) ## gives the levels, i.e. 'TRUE' means daylight, 'FALSE' means nighttime and darkness
## [1] FALSE  TRUE

Check for temporal outliers

Check if there are strange dates, or dates before you collared the animal (e.g., usually the collar is activated and tested before the animal is collared):

table(dat_anim_na$date)
## 
## 2018-10-30 2018-10-31 2018-11-01 2018-11-02 2018-11-03 2018-11-04 2018-11-05 
##         24         40         35         39         42         34         34 
## 2018-11-06 2018-11-07 2018-11-08 2018-11-09 2018-11-10 2018-11-11 2018-11-12 
##         41         42         42         41         38         38         39 
## 2018-11-13 2018-11-14 2018-11-15 2018-11-16 2018-11-17 2018-11-18 2018-11-19 
##         41         38         34         38         44         34         35 
## 2018-11-20 2018-11-21 2018-11-22 2018-11-23 2018-11-24 2018-11-25 2018-11-26 
##         36         34         40         46         40         37         41 
## 2018-11-27 2018-11-28 2018-11-29 2018-11-30 2018-12-01 2018-12-02 2018-12-03 
##         43         39         34         42         45         37         39 
## 2018-12-04 2018-12-05 2018-12-06 2018-12-07 2018-12-08 2018-12-09 2018-12-10 
##         38         39         38         31         38         36         36 
## 2018-12-11 2018-12-12 2018-12-13 2018-12-14 2018-12-15 2018-12-16 2018-12-17 
##         39         34         38         37         37         43         43 
## 2018-12-18 2018-12-19 2018-12-20 2018-12-21 2018-12-22 2018-12-23 2018-12-24 
##         37         31         41         36         38         42         44 
## 2018-12-25 2018-12-26 2018-12-27 2018-12-28 2018-12-29 2018-12-30 2018-12-31 
##         32         30         38         38         36         37         39 
## 2019-01-01 2019-01-02 2019-01-03 2019-01-04 2019-01-05 2019-01-06 2019-01-07 
##         33         43         37         36         39         39         45 
## 2019-01-08 2019-01-09 2019-01-10 2019-01-11 2019-01-12 2019-01-13 2019-01-14 
##         36         38         41         33         26         35         40 
## 2019-01-15 2019-01-16 2019-01-17 2019-01-18 2019-01-19 2019-01-20 2019-01-21 
##         28         35         35         40         39         34         31 
## 2019-01-22 2019-01-23 2019-01-24 2019-01-25 2019-01-26 2019-01-27 2019-01-28 
##         36         36         34         36         43         34         33 
## 2019-01-29 2019-01-30 2019-01-31 2025-12-26 
##         35         34         27          1

It might be easier to spot visually. We use the {ggplot2} package here as it recognizes and respects the date format:

ggplot(dat_anim_na, aes(date)) +
  geom_bar() +
  theme_bw()

## compare with plot(table(dat_anim_na$date))

There is a strange date → 2025-12-26!

Delete this data row:

delme <- which(dat_anim_na$date == "2025-12-26")
dat_anim_na[delme,]                 ## check observation
##      key.bin.checksum tag.serial.number     start.timestamp longitude latitude
## 2622       2948730959              5334 2025-12-26 04:00:00  13.47851 52.49133
##      height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 2622                   78.8           3      A                   23
##             timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 2622 2025-12-26 04:00:24.000            3582                3287           3
##      speed.over.ground heading.degree speed.accuracy.estimate
## 2622              0.03              0                    0.64
##      horizontal.accuracy.estimate yearday month hour kweek       date  sunrise
## 2622                         4.35     360    12    4    52 2025-12-26 8.390617
##        sunset daylength daytime
## 2622 15.80777  7.417151   FALSE
dat_anim_na <- dat_anim_na[-delme,] ## delete the strange date and 
table(dat_anim_na$date)             ## check again
## 
## 2018-10-30 2018-10-31 2018-11-01 2018-11-02 2018-11-03 2018-11-04 2018-11-05 
##         24         40         35         39         42         34         34 
## 2018-11-06 2018-11-07 2018-11-08 2018-11-09 2018-11-10 2018-11-11 2018-11-12 
##         41         42         42         41         38         38         39 
## 2018-11-13 2018-11-14 2018-11-15 2018-11-16 2018-11-17 2018-11-18 2018-11-19 
##         41         38         34         38         44         34         35 
## 2018-11-20 2018-11-21 2018-11-22 2018-11-23 2018-11-24 2018-11-25 2018-11-26 
##         36         34         40         46         40         37         41 
## 2018-11-27 2018-11-28 2018-11-29 2018-11-30 2018-12-01 2018-12-02 2018-12-03 
##         43         39         34         42         45         37         39 
## 2018-12-04 2018-12-05 2018-12-06 2018-12-07 2018-12-08 2018-12-09 2018-12-10 
##         38         39         38         31         38         36         36 
## 2018-12-11 2018-12-12 2018-12-13 2018-12-14 2018-12-15 2018-12-16 2018-12-17 
##         39         34         38         37         37         43         43 
## 2018-12-18 2018-12-19 2018-12-20 2018-12-21 2018-12-22 2018-12-23 2018-12-24 
##         37         31         41         36         38         42         44 
## 2018-12-25 2018-12-26 2018-12-27 2018-12-28 2018-12-29 2018-12-30 2018-12-31 
##         32         30         38         38         36         37         39 
## 2019-01-01 2019-01-02 2019-01-03 2019-01-04 2019-01-05 2019-01-06 2019-01-07 
##         33         43         37         36         39         39         45 
## 2019-01-08 2019-01-09 2019-01-10 2019-01-11 2019-01-12 2019-01-13 2019-01-14 
##         36         38         41         33         26         35         40 
## 2019-01-15 2019-01-16 2019-01-17 2019-01-18 2019-01-19 2019-01-20 2019-01-21 
##         28         35         35         40         39         34         31 
## 2019-01-22 2019-01-23 2019-01-24 2019-01-25 2019-01-26 2019-01-27 2019-01-28 
##         36         36         34         36         43         34         33 
## 2019-01-29 2019-01-30 2019-01-31 
##         35         34         27
plot(table(dat_anim_na$date))       ## plot the number of fixes per day

Save the cleaned file

Finally, we save the processed data file that we will use for exploration and analysis into the subfolder .../output/data-proc. There are two options we can use:

  • R data file - can only be opened/ read with R by using function readRDS()
saveRDS(dat_anim_na, file = here("output",  "tag5334_gps_proc.Rds"))
  • interchange file format .csv
write.csv(dat_anim_na, file = here("output",  "tag5334_gps_proc.csv"))

Check your output folder for these files. The efficient ‘.Rds’ file storage is about 3 times smaller than the ‘.csv’ file.

Exercise 4.1

This is a recap from Course2. Please plot the animal relocations in two different colours based on the column daytime, using one of the options to create maps with {ggplot2}, {ggmap}, {leaflet} , {tmap} or {ggspatial}. Use the processed file of fox Q (tag5334_gps_proc) and do it in a separate script. Save your script as Course4_Exercise1_*yourname*.R.

Hint: Remember to load the relevant libraries. Note that you might want to plot the locations in the correct spatial dimensions by projecting it using the functions st_as_sf() and st_transform().



Data exploration

Load the cleaned data file

I recommend that you store the raw file safely and continue with the cleaned and processed file after major data manipulations have been conducted to minimize errors. I’d even suggest to make these steps in different R-scripts.

anim_proc <- readRDS(file = here("output", "tag5334_gps_proc.Rds"))
# head(anim_proc)

Have a look whether there are gaps/ missing days in the data, or whether we have approximately a regular number of fixes each day. We can use {ggplot2} to plot true dates:

ggplot(anim_proc, aes(date)) +
  geom_bar() +
  theme_bw()

# mind: is the same as already plotted above:
# plot(table(anim_proc$date))

Now let’s plot the Julian day, yearday:

plot(table(anim_proc$yearday))  ## what is the difference to the plot above?

The number of fixes (locations taken) per day seems to be quite regular. Mind the difference when plotting per date (sorted) and yearday (1-365), irrespective of the year.

Plot activity

Now we can have a look at the distribution of the fixes during the day, which is a hint on the active phase of the animal:

plot(table(anim_proc$hour))   ## plot the number of fixes per hour

More fixes during the night - this is the active phase. As these are circular data, we can plot the number of fixes as sign for activity (knowing that our collars did not record when foxes were inactive! Otherwise, we cannot distinguish missing data from inactive times! -> bias).

## simple circular plot (rose diagram)
timetoplot <- circular(anim_proc$hour %% 24, ## convert to 24 hrs = bins
                       units = "hours", template = "clock24")

## Note: use namespace `circular::` here as there are multiple functions called `rose.diag()`
circular::rose.diag(timetoplot, bin = 24, col = "blue",
                    main = "Events by Hour (sqrt scale)", prop = 3)

Check whether the activity during the day was before sunrise (our ‘daytime’ column, which is either TRUE (=yes, daylight) or FALSE (dark hours):

# with ggplot
# code adapted from https://gist.github.com/mattbaggott/4361381
ggplot(anim_proc, aes(x = hour,fill = daytime)) +
  geom_bar(width = 1, color = "grey20") +
  coord_polar(start = -0.15) + ## to center setchange start
  theme_minimal() +
  scale_x_continuous(breaks = 0:23) +
  scale_fill_brewer(palette = "Set2", name = NULL, labels = c("Night", "Day")) +
  labs(x = NULL, y = "Count", title = "Events by Time of the Day")

ggsave(filename = 'plot_activity.png', path=here('plots'), 
       width = 10, height = 10, units = "cm")

Apart from the outliers between 11 and 13 o’clock, the active phase was during the dark hours.

Actogramm

Convert the data into a spatial object

Now we transform the data into sf and sp objects and assign the coordinate reference system:

# transform into spatial simple feature sf object
mydf_sf <- st_as_sf(x = data.frame(anim_proc),
                       coords = c("longitude", "latitude"),
                       crs = 4326,
                       sf_column_name = "geometry" )

# transform into SpatialPointsDataFrame  - for crosschecking
mydf_sp <- as(mydf_sf, "Spatial")

And then we project the reference system from angular units to a planar coordinate reference system in meters:

# transform CRS to projected one in meter distance units
mydf_sf_trans <-  st_transform(mydf_sf, 3035 )  # EPSG-code  
mydf_sp_trans <-  spTransform(mydf_sp, CRS("+init=epsg:3035")) ## warning here, as 'sp' gets outdated

Recently, there are issues with the missing datum in the CRS-specifications. We ignore this for now. More on the issue of moving from proj4 to proj6 in the future:
* https://inbo.github.io/tutorials/tutorials/spatial_crs_coding/
* https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html

Make a quick plot of the data

Recap for styling the plot:
https://www.r-bloggers.com/2021/12/introduction-to-geospatial-visualization-with-the-tmap-package/

tmap_mode(mode = "view")

tm_shape(shp = mydf_sf_trans) + 
  tm_dots(size = 0.01, 
          col = "daytime",  
          alpha = 0.5)

Have a deep look at the data: did the fox really swim? (points in Lake Rummelsburg). Or are these outliers? No: the lake was frozen - check the date! Note: you really need to know your animals and the area before analysing data!

Basic movement metrics

Transform your data into an {amt} object

Check here for the possibilities of what to do with the {amt}-package:
* https://cran.r-project.org/web/packages/amt/vignettes/p1_getting_started.html

Because we are dealing with spatial measures (e.g. area, distance, perimeter,…), I recommend to work with projected coordinates. To create a move-object, the coordinates of the projected CRS need to be appended to a data.frame:

## assign columns with coordinates_epsg to add later into the data frame for amt-object
mydf_sf_trans$X_3035 <- st_coordinates(mydf_sf_trans)[,1]
mydf_sf_trans$Y_3035 <- st_coordinates(mydf_sf_trans)[,2]

## new object as data frame
anim_for_move <- st_drop_geometry(mydf_sf_trans)
head(anim_for_move)
##   key.bin.checksum tag.serial.number     start.timestamp height.above.ellipsoid
## 1       2135768252              5334 2018-10-30 04:00:00                   76.3
## 2       2210577241              5334 2018-10-30 04:20:00                   81.6
## 3       1824790132              5334 2018-10-30 04:40:00                   66.5
## 4       4079427513              5334 2018-10-30 07:00:00                   73.5
## 5        380027806              5334 2018-10-30 08:00:00                   71.1
## 6       3352806116              5334 2018-10-30 12:00:00                   71.5
##   type.of.fix status used.time.to.get.fix        timestamp.of.fix
## 1           3      A                   22 2018-10-30 04:00:23.000
## 2           3      A                   18 2018-10-30 04:20:19.000
## 3           3      A                   53 2018-10-30 04:40:54.000
## 4           3      A                   24 2018-10-30 07:00:25.000
## 5           3      A                   47 2018-10-30 08:00:48.000
## 6           3      A                   23 2018-10-30 12:00:24.000
##   battery.voltage fix.battery.voltage temperature speed.over.ground
## 1            3604                3389          10              1.64
## 2            3616                3333          11              4.43
## 3            3624                3374          13              0.07
## 4            3659                3491          11              0.01
## 5            3660                3470          13              0.04
## 6            3664                3508          14              0.02
##   heading.degree speed.accuracy.estimate horizontal.accuracy.estimate yearday
## 1         233.31                    1.67                         7.94     303
## 2         327.02                    3.79                        59.65     303
## 3           0.00                    0.73                        15.36     303
## 4           0.00                    0.65                         7.17     303
## 5           0.00                    0.89                         8.45     303
## 6           0.00                    0.44                         3.84     303
##   month hour kweek       date  sunrise   sunset daylength daytime  X_3035
## 1    10    4    44 2018-10-30 7.033396 16.62727  9.593878   FALSE 4556608
## 2    10    4    44 2018-10-30 7.033335 16.62719  9.593855   FALSE 4556681
## 3    10    4    44 2018-10-30 7.033371 16.62723  9.593855   FALSE 4556645
## 4    10    7    44 2018-10-30 7.033367 16.62722  9.593855   FALSE 4556649
## 5    10    8    44 2018-10-30 7.033370 16.62722  9.593853    TRUE 4556646
## 6    10   12    44 2018-10-30 7.033371 16.62722  9.593852    TRUE 4556646
##    Y_3035
## 1 3270698
## 2 3270731
## 3 3270729
## 4 3270729
## 5 3270731
## 6 3270732
## now turn into amt object

all(complete.cases(anim_for_move))
## [1] TRUE
my_fox <- make_track(tbl = anim_for_move, 
                     .x  = X_3035, 
                     .y  = Y_3035, 
                     .t  = start.timestamp, 
                     id = tag.serial.number, ## has to be a column of data set
                     crs = 3035)             ## has to be set as EPSG code
                    # all_cols = TRUE)        ## if all columns need to be kept

head(my_fox)
## # A tibble: 6 × 4
##         x_       y_ t_                     id
## *    <dbl>    <dbl> <dttm>              <int>
## 1 4556608. 3270698. 2018-10-30 04:00:00  5334
## 2 4556681. 3270731. 2018-10-30 04:20:00  5334
## 3 4556645. 3270729. 2018-10-30 04:40:00  5334
## 4 4556649. 3270729. 2018-10-30 07:00:00  5334
## 5 4556646. 3270731. 2018-10-30 08:00:00  5334
## 6 4556646. 3270732. 2018-10-30 12:00:00  5334
## check for duplicated time stamps
any(duplicated(my_fox$t_))
## [1] FALSE
## FALSE - great, let's go!

A hint (for the exercises later, so please read!): check the argument ‘id’. Here, I put one tag-id number, however, if multiple individuals are used, you can provide a column with multiple IDs.

Movement metrics

## Time lags between the fixes:
summarize_sampling_rate(my_fox) ## varies a lot, but median here 20 min, max 4 hrs
## # A tibble: 1 × 9
##     min    q1 median  mean    q3   max    sd     n unit 
##   <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1  17.2    20     20  38.5    20   240  46.8  3502 min
## The distance ranges form 0.5 m to 1 km,
## but we have here the 4 hrs time lags included
range(step_lengths(my_fox), na.rm =TRUE)
## [1]    0.5468472 1028.7766731
## speed
class(my_fox)
## [1] "track_xyt"  "track_xy"   "tbl_df"     "tbl"        "data.frame"
#max(amt::speed(my_fox),na.rm=TRUE) ## 0.8, measured in m/s

## Let's look only at 60 min intervals:
my_fox_60 <- track_resample(my_fox, rate = minutes(60), tolerance = minutes(5))
range(step_lengths(my_fox_60), na.rm =TRUE)
## [1]    0.7395256 1184.1579983
## can you explain why the step length gets larger?
## what happens with speed?
#max(speed(my_fox_60),na.rm=TRUE)

Turning angles:

The calculation of turning angles only makes sense with very high resolution data. With 20 min intervals between fixes, turning angles do not make much sense any more, as the animals can have turned and moved a lot in between. I nevertheless show the code, because turning angles are often used to describe behavior.
Signer et al. 2018; DOI: 10.1002/ece3.4823
We will also choose to keep only those bursts (subsets of the track with constant sampling rate, within the specified tolerance) with at least three relocations, the minimum required to calculate a turn angle (amt::filter_min _ n _ burst). The following code implements those choices and translates from a point representation to a step (step length, turn angle) representation of the data.

## first transform track into steps
tmp1 <- amt::filter_min_n_burst(my_fox_60, min_n=3)

## then derive turning angles - column 'ta_' will be appended
fox_steps_60 <- steps_by_burst(tmp1) ## turn. ang. from -pi to pi = -180° to 180°
range(as_degree(fox_steps_60$ta_), na.rm = TRUE) ## convert to degree
## [1] -179.9449  179.9833
hist(as_degree(fox_steps_60$ta_)) ## as said, with hourly interval, this is not exciting

The homerange concept - MCP, kernelUD (aKDE)

See lecture!

Minimum Convex Polygon

Get the home range size.
https://cran.r-project.org/web/packages/amt/vignettes/p2_hr.html

mcp_fox_95 <- amt::hr_mcp(x = my_fox, ## the track object
                        levels = c(0.95)) # define the level of the mcp. 
plot(mcp_fox_95) ## quick and uggly

## if you are not sure about how to get the info of the object, use str()
mcp_fox_95$mcp$area/ 1e6 ## area converted to km2
## 0.586523 [m^2]
## same as 
hr_area(mcp_fox_95)
## # A tibble: 1 × 3
##   level what        area
##   <dbl> <chr>      <dbl>
## 1  0.95 estimate 586523.
## calculated MCP at different levels
mcp_fox <- hr_mcp(my_fox, levels = seq(0.2, 1, 0.1))
hr_area(mcp_fox)
## # A tibble: 9 × 3
##   level what        area
##   <dbl> <chr>      <dbl>
## 1   1   estimate 684232.
## 2   0.9 estimate 509770.
## 3   0.8 estimate 305613.
## 4   0.7 estimate 232182.
## 5   0.6 estimate 162707.
## 6   0.5 estimate 131338.
## 7   0.4 estimate 123252.
## 8   0.3 estimate 110913.
## 9   0.2 estimate  79550.
plot(hr_area(mcp_fox)$area ~ hr_area(mcp_fox)$level, type = 'l')

Convert to sf object to store as shapefile, for example, and make a plot 

mcp_fox_95_sf <- sf::st_as_sf(mcp_fox_95$mcp)

st_write(mcp_fox_95_sf,  
         dsn = here("output", "mcp_fox_95_sf.shp"),
         delete_layer = TRUE)
## Deleting layer `mcp_fox_95_sf' using driver `ESRI Shapefile'
## Writing layer `mcp_fox_95_sf' to data source 
##   `M:/IZW/github_local_Mdrive/d6_teaching_collection/output/mcp_fox_95_sf.shp' using driver `ESRI Shapefile'
## Writing 1 features with 3 fields and geometry type Polygon.
ggplot() + 
  geom_sf(data = mcp_fox_95_sf) + 
  geom_sf(data = mydf_sf_trans, size=0.1, col = 'red') + ## df converted to sf above
  coord_sf()

Kernel utility density

Calculate the density kernel:

kde_fox_95 <- amt::hr_kde(x = my_fox, levels = c(0.95))
plot(kde_fox_95)

hr_area(kde_fox_95)
## # A tibble: 1 × 3
##   level what        area
##   <dbl> <chr>      <dbl>
## 1  0.95 estimate 532654.
hr_area(mcp_fox_95) ## for comparison
## # A tibble: 1 × 3
##   level what        area
##   <dbl> <chr>      <dbl>
## 1  0.95 estimate 586523.
## spatial density distribution
# plot(kde_fox_95$ud,col =viridis(100, direction = -1)) ## density distribution
# contour(kde_fox_95$ud)

Also here, convert to sf object to store as shapefile, for example, and make a plot 

kde_fox_95_sf <- hr_isopleths(kde_fox_95) ## nicely, this is already an sf object

st_write(kde_fox_95_sf, 
         dsn = here("output", "kde_fox_95_sf.shp"),
         delete_layer = TRUE)
## Deleting layer `kde_fox_95_sf' using driver `ESRI Shapefile'
## Writing layer `kde_fox_95_sf' to data source 
##   `M:/IZW/github_local_Mdrive/d6_teaching_collection/output/kde_fox_95_sf.shp' using driver `ESRI Shapefile'
## Writing 1 features with 3 fields and geometry type Multi Polygon.
ggplot() + 
  geom_sf(data = mcp_fox_95_sf, alpha=0.5) + 
  geom_sf(data = kde_fox_95_sf, alpha=0.1) + 
  geom_sf(data = mydf_sf_trans, size=0.1, col = 'red') + ## df converted to sf above
  coord_sf()

Resource selection

Nothing makes sense but in the light of environmental information…

https://cran.r-project.org/web/packages/amt/vignettes/p3_rsf.html

Let’s first thin and filter the dataset to one location per day

## Let's look only at 60 min intervals:
my_fox_24 <- track_resample(my_fox, rate = hours(24), tolerance = minutes(25))
dim(my_fox_24)
## [1] 95  5

Environmental covariates

First, load some environmental data. Check what is available in your data_berlin-course folder, check for the CRS and assign if missing:

maps_wd   <- here("data","data_berlin","geo_raster_current_gtif")  
lf <- list.files(maps_wd, full.names = TRUE); lf
## [1] "M:/IZW/github_local_Mdrive/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/imperviousness_2012_100m_3035.tif"    
## [2] "M:/IZW/github_local_Mdrive/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/noise_daynight_2017_100m_3035.tif"    
## [3] "M:/IZW/github_local_Mdrive/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/summer_temp_100m_3035.tif"            
## [4] "M:/IZW/github_local_Mdrive/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/tree_cover_density_2012_100m_3035.tif"
## [5] "M:/IZW/github_local_Mdrive/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/water_bodies_2010_100m_3035.tif"
b_imperv <- terra::rast(lf[1]) 
b_sutemp <- terra::rast(lf[3]) 
b_forest <- terra::rast(lf[4])
b_noise  <- terra::rast(lf[2])

water_raster  <- rast(here::here("data","data_berlin","geo_raster_current_gtif","water_bodies_2010_100m_3035.tif"))

## chech CRS:
b_sutemp ## coord. ref. is assigned because it is a geoTiff, not ascii-file.
## class       : SpatRaster 
## dimensions  : 570, 657, 1  (nrow, ncol, nlyr)
## resolution  : 100, 100  (x, y)
## extent      : 4521040, 4586740, 3243800, 3300800  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source      : summer_temp_100m_3035.tif 
## name        : summer_temp_100m_3035 
## min value   :                  2052 
## max value   :                  3133
terra::plot(b_sutemp)

## combine to raster stack
## the variables are on very different scales - we will scale them to have standardized values
b_env_var_stack <- terra::rast(list(scale(b_imperv),
                                    scale(b_sutemp), 
                                    scale(b_forest),
                                    scale(b_noise)))
names(b_env_var_stack) <- c("imp","stmp","tree","noise")

terra::plot(b_env_var_stack) 

## single plot with points
#terra::plot(b_env_var_stack[[4]])
#points(my_fox_24, col = 'red')

Compute a correlation matrix before you enter all variables into the model to check for multicollinearity

M <- terra::layerCor(b_env_var_stack, fun = "pearson")
corrplot(M$correlation, type="upper",tl.col="black", tl.srt=45)

Now create random points in MCP home range of the fox. Default is 10 times more random points than points. It automatically sets a column ‘case_’ = FALSE, because it is a random point, not a true animal location.

rp_fox_24 <- random_points(my_fox_24, n=dim(my_fox_24)[1])
plot(rp_fox_24)

Now we extract the covariates (several ways to do so, see Course 2 Spatial R and function ‘st_intersection()’ or ‘extract’, but we use the function ‘extract_covariates()’ of the amt-package here)

env_rp_fox_24 <- amt::extract_covariates(x=rp_fox_24, covariates=b_env_var_stack)
head(env_rp_fox_24)
## # A tibble: 6 × 7
##   case_       x_       y_      imp  stmp   tree   noise
## * <lgl>    <dbl>    <dbl>    <dbl> <dbl>  <dbl>   <dbl>
## 1 FALSE 4557255. 3270205. -0.00158 -1.59  0.940 -0.264 
## 2 FALSE 4556834. 3270444.  0.0240  -1.25 -0.197 -0.0715
## 3 FALSE 4557119. 3270413. -0.219   -1.51  0.403 -0.336 
## 4 FALSE 4557407. 3270298.  0.942   -1.69  0.965 -0.462 
## 5 FALSE 4557087. 3270622.  1.15    -1.61 -2.19  -0.264 
## 6 FALSE 4557308. 3270349.  0.544   -1.67  1.05  -0.325

We do the same with the fox locations. Note that we have a different structure of the data.frame with more columns. And we have to add the column ‘case_’

env_my_fox_24 <- extract_covariates(x=my_fox_24, covariates=b_env_var_stack)
head(env_my_fox_24)
## # A tibble: 6 × 9
##         x_       y_ t_                     id burst_    imp  stmp    tree  noise
## *    <dbl>    <dbl> <dttm>              <int>  <dbl>  <dbl> <dbl>   <dbl>  <dbl>
## 1 4556608. 3270698. 2018-10-30 04:00:00  5334      1  1.18  -1.37 -0.721  -0.148
## 2 4556644. 3270734. 2018-10-31 03:40:00  5334      1  1.01  -1.54 -0.0547 -0.160
## 3 4556808. 3270550. 2018-11-01 03:40:00  5334      1  0.561 -1.36  0.0990 -0.443
## 4 4557182. 3270347. 2018-11-02 03:20:00  5334      1 -0.865 -1.57  0.734  -0.213
## 5 4557442. 3270296. 2018-11-03 03:00:00  5334      1  1.36  -1.80 -2.19   -0.438
## 6 4556524. 3270576. 2018-11-04 02:40:00  5334      1  1.04  -1.15 -2.19    0.277
env_my_fox_24$case_ <- 'TRUE'
names(env_my_fox_24) ## select same columns as random point df
##  [1] "x_"     "y_"     "t_"     "id"     "burst_" "imp"    "stmp"   "tree"  
##  [9] "noise"  "case_"
env_my_fox_24 <- env_my_fox_24[,c(10,1,2,6:9)]
head(env_my_fox_24)  
## # A tibble: 6 × 7
##   case_       x_       y_    imp  stmp    tree  noise
## * <chr>    <dbl>    <dbl>  <dbl> <dbl>   <dbl>  <dbl>
## 1 TRUE  4556608. 3270698.  1.18  -1.37 -0.721  -0.148
## 2 TRUE  4556644. 3270734.  1.01  -1.54 -0.0547 -0.160
## 3 TRUE  4556808. 3270550.  0.561 -1.36  0.0990 -0.443
## 4 TRUE  4557182. 3270347. -0.865 -1.57  0.734  -0.213
## 5 TRUE  4557442. 3270296.  1.36  -1.80 -2.19   -0.438
## 6 TRUE  4556524. 3270576.  1.04  -1.15 -2.19    0.277

Combine datasets for analysis and convert binary TRUE/FALSE (column ‘case_’) into numeric.

df_rsf_fox_24 <- rbind(env_rp_fox_24,env_my_fox_24)
head(df_rsf_fox_24)
## # A tibble: 6 × 7
##   case_       x_       y_      imp  stmp   tree   noise
## * <chr>    <dbl>    <dbl>    <dbl> <dbl>  <dbl>   <dbl>
## 1 FALSE 4557255. 3270205. -0.00158 -1.59  0.940 -0.264 
## 2 FALSE 4556834. 3270444.  0.0240  -1.25 -0.197 -0.0715
## 3 FALSE 4557119. 3270413. -0.219   -1.51  0.403 -0.336 
## 4 FALSE 4557407. 3270298.  0.942   -1.69  0.965 -0.462 
## 5 FALSE 4557087. 3270622.  1.15    -1.61 -2.19  -0.264 
## 6 FALSE 4557308. 3270349.  0.544   -1.67  1.05  -0.325
df_rsf_fox_24$case_num <- ifelse(df_rsf_fox_24$case_ == 'FALSE',0,1)

RSF model fit

Scaling variables, so that they all have the same range of values, can be helpful when plotting - only few plot functions below have inbuilt functions to standardize while plotting.

rsf_fox <- glm(data    = df_rsf_fox_24, 
               #formula = case_num ~ scale(noise) + scale(imp) + scale(stmp) + scale(tree), 
               formula = case_num ~ noise + imp + stmp + tree, 
               family  = binomial(link = "logit"))
summary(rsf_fox)
## 
## Call:
## glm(formula = case_num ~ noise + imp + stmp + tree, family = binomial(link = "logit"), 
##     data = df_rsf_fox_24)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2651     1.7478  -0.152 0.879451    
## noise         3.7353     1.0207   3.659 0.000253 ***
## imp           1.4602     0.2876   5.078 3.82e-07 ***
## stmp         -0.7634     1.2123  -0.630 0.528877    
## tree          1.0566     0.2421   4.363 1.28e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 362.81  on 284  degrees of freedom
## Residual deviance: 305.37  on 280  degrees of freedom
## AIC: 315.37
## 
## Number of Fisher Scoring iterations: 4
## model frame - if you need to check 
#model.frame(rsf_fox)
#nrow(model.frame(rsf_fox)) == nrow(rsf_fox)

Effects plots

Plot the estimates as effects plots and forest plots.

plot(effects::allEffects(rsf_fox),lty=1)

Forest plot

https://www.jtools.jacob-long.com/articles/summ.html

plot_summs(rsf_fox, inner_ci_level = .75) #, plot.distributions = TRUE

The r-package `{flextable}´ provides an easy way to export the model summary from R into an excel or word file format.
To sum up: The results - that a fox prefers houses and noise - do not really make sense. What could be the reason be?
You could repeat the analysis with day and nighttime locations of the fox….


Predict the model

The predict function will calculate the response on the logit scale, therefore we have to back-transform the values with the inverse logit to a probability.

f_prob  <- function(var2) {1/ (1 + exp(x = -var2))}  # inverse logit

## the prediction is made on logit scale
fox_pred_logit <- terra::predict(object = b_env_var_stack, 
                                model = rsf_fox)

## transform logit to probability (use inverse logit function above)
fox_pred_prob <- f_prob(fox_pred_logit)
# plot map
terra::plot(fox_pred_prob, col = viridis(100) )
terra::plot(water_raster, col = "darkslategray1",  legend=FALSE, add = TRUE)

# writeRaster(fox_pred, here("outpu"))

-> go to Exercise 4.2;
check argument ‘nest: -id’ in amt-package to work with several animals or make a subset for ach animal

Coming soon: Recurse analysis, HMM, …

END

Adding the session info can be very helpful when going back to old scripts or using scripts of others:

Session Info 
sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8   
## [3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=C                      
## 
## time zone: Europe/Berlin
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] solartime_0.0.4   corrplot_0.95     viridis_0.6.5     viridisLite_0.4.2
##  [5] jtools_2.3.0      effects_4.2-2     carData_3.0-5     ctmm_1.2.0       
##  [9] amt_0.2.2.0       move_4.2.6        raster_3.6-31     sp_2.2-0         
## [13] geosphere_1.5-20  circular_0.5-1    tmap_4.0          ggplot2_3.5.1    
## [17] terra_1.8-21      sf_1.0-19         lubridate_1.9.4   here_1.0.1       
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3      rstudioapi_0.17.1       jsonlite_1.8.9         
##   [4] wk_0.9.4                magrittr_2.0.3          estimability_1.5.1     
##   [7] farver_2.1.2            nloptr_2.1.1            rmarkdown_2.29         
##  [10] ragg_1.3.3              vctrs_0.6.5             memoise_2.0.1          
##  [13] minqa_1.2.8             base64enc_0.1-3         htmltools_0.5.8.1      
##  [16] forcats_1.0.0           leafsync_0.1.0          survey_4.4-2           
##  [19] broom_1.0.7             s2_1.1.7                sass_0.4.9             
##  [22] parallelly_1.40.1       KernSmooth_2.23-24      bslib_0.8.0            
##  [25] htmlwidgets_1.6.4       stars_0.6-7             cachem_1.1.0           
##  [28] lifecycle_1.0.4         pkgconfig_2.0.3         cols4all_0.8           
##  [31] Matrix_1.7-1            R6_2.6.1                fastmap_1.2.0          
##  [34] rbibutils_2.3           future_1.34.0           digest_0.6.37          
##  [37] colorspace_2.1-1        furrr_0.3.1             rprojroot_2.0.4        
##  [40] leafem_0.2.3            textshaping_0.4.1       crosstalk_1.2.1        
##  [43] lwgeom_0.2-14           labeling_0.4.3          spacesXYZ_1.5-1        
##  [46] timechange_0.3.0        httr_1.4.7              abind_1.4-8            
##  [49] compiler_4.4.2          proxy_0.4-27            withr_3.0.2            
##  [52] pander_0.6.5            backports_1.5.0         DBI_1.2.3              
##  [55] logger_0.4.0            broom.mixed_0.2.9.6     MASS_7.3-61            
##  [58] tmaptools_3.1-1         leaflet_2.2.2           classInt_0.4-11        
##  [61] tools_4.4.2             units_0.8-5             leaflegend_1.2.1       
##  [64] nnet_7.3-19             glue_1.8.0              leafgl_0.2.2           
##  [67] nlme_3.1-166            grid_4.4.2              checkmate_2.3.2        
##  [70] generics_0.1.3          gtable_0.3.6            leaflet.providers_2.0.0
##  [73] class_7.3-22            rmdformats_1.0.4        tidyr_1.3.1            
##  [76] data.table_1.16.4       utf8_1.2.4              xml2_1.3.6             
##  [79] pillar_1.10.1           mitools_2.4             splines_4.4.2          
##  [82] dplyr_1.1.4             lattice_0.22-6          survival_3.7-0         
##  [85] tidyselect_1.2.1        knitr_1.49              jsonify_1.2.2          
##  [88] gridExtra_2.3           bookdown_0.41           xfun_0.49              
##  [91] yaml_2.3.10             boot_1.3-31             evaluate_1.0.1         
##  [94] codetools_0.2-20        tibble_3.2.1            cli_3.6.3              
##  [97] systemfonts_1.1.0       Rdpack_2.6.2            munsell_0.5.1          
## [100] jquerylib_0.1.4         dichromat_2.0-0.1       Rcpp_1.0.14            
## [103] globals_0.16.3          png_0.1-8               XML_3.99-0.17          
## [106] parallel_4.4.2          lme4_1.1-35.5           listenv_0.9.1          
## [109] mvtnorm_1.3-3           scales_1.3.0            e1071_1.7-16           
## [112] insight_1.0.2           purrr_1.0.2             rlang_1.1.4